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Statement  of  Problem  Studied 


Computational  Model  for  Domain  Structure  Evolution  in  Ferroelectrics 

Final  Report 
PI:  Chad  M.  Landis 

Students:  Wenyuan  Li  and  Dorinamaria  Carka 
Post-doc:  Antonios  Kontsos 

The  goal  of  the  proposed  research  is  to  develop  a  sound  theoretical  framework  and 
advanced  numerical  techniques  to  simulate  the  evolution  of  domain  microstructures  in 
ferroelectric  materials.  Ferroelectric  ceramics  are  presently  being  used  in  a  broad  range 
of  applications  including  fuel  injectors  for  high  efficiency- low  emission  diesel  engines, 
actuators  for  active  control  of  helicopter  rotor  blades  and  underwater  vehicle  control 
surfaces,  and  ultrasonic  rotary  inchworm  motors  with  high  power  and  torque  densities. 
Additionally,  ferroelectric  thin  films  are  used  for  data  storage  in  non-volatile 
ferroelectric  random  access  memory  (NVFRAM),  sensing  and  actuation  in 
microelectromechanical  system  (MEMS),  and  in  nonlinear  optics.  An  understanding  of 
microstructural  evolution  and  domain  dynamics  is  necessary  for  further  development  of 
micro/nano- ferroelectric  device  technology.  We  developed  a  combined  theoretical  and 
numerical  modeling  framework  in  order  to  investigate  the  interactions  of  domain  walls 
with  surfaces,  grain  boundaries,  charges,  dislocations  and  other  types  of  defects.  Such 
studies  will  aid  in  the  understanding  of  device  performance  and  of  material  failure 
mechanisms. 


Summary  of  Important  Results 


Theoretical  Foundations 

Here  we  review  the  fundamental  governing  phase-field  equations.  Standard  index 
notation  is  used  with  summation  implied  over  repeated  indices,  the  overdot  represents 
the  first  derivative  with  respect  to  time,  and  ,j  represents  partial  differentiation  with 
respect  to  the  x  coordinate  direction. 

The  static  equilibrium  equations  in  any  arbitrary  volume  V  and  its  bounding  surfaces  S 


yield, 

<7 .  +b 

Jl,J  * 

=  0  in  V 

(1) 

a..  =  a 

ij  ji 

in  V 

(2) 

g  n.  =  t 

v  j 

on  S 

l 

(3) 

Where  ov  are  the  Cartesian  components  of  the  Cauchy  stress,  b.  the  components  of  a 
body  force  per  unit  volume,  n  the  components  of  a  unit  vector  normal  to  a  surface 
element,  and  t  the  tractions  applied  to  a  surface.  Under  the  assumptions  of  linear 
kinematics,  the  strain  components  £  are  related  to  the  displacements  u.  as 


£.. 
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1/ 
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[u.  +  u.  din  V 

V  hj  j,t/ 


(4) 


The  electrical  quantities  of  electric  field,  E.  ,  electric  potential,  4> ,  electric  displacement, 
D  ,  volume  charge  density,  q ,  and  surface  charge  density,  u  ,  are  governed  by  the 
quasi-static  forms  of  Maxwell’s  equations.  Specifically, 

D  .  —  q  =  0  in  V  (5) 

Dn  =  —cu  on  S  (6) 

E,  =  -t  (7) 

within  the  phase-field  modeling  approach  the  free  energy  will  also  be  required  to  depend 
on  an  order  parameter  and  its  gradient.  Mathematically,  the  order  parameter  is  used  to 
describe  the  different  variant  types,  i.e.  the  spontaneous  states  that  are  possible.  For 
the  case  of  ferroelectrics,  the  physically  natural  order  parameter  is  the  electrical 


polarization  P .  Note  that  the  relationship  between  electric  field,  electric  displacement 
and  polarization  in  matter  is  given  as 


D  =  P  +  k,,E 


(8) 


Here  kq  is  the  permittivity  of  free  space.  Given  that  the  free  energy  is  permitted  to 

depend  on  an  order  parameter,  a  set  of  micro- forces  are  introduced  that  are  work- 
conjugate  to  the  order  parameter  and  its  gradient.  The  micro- force  balance  is  written 
as, 


£  .  +  7 r.  +  7.  =  0 


(9) 


where  we  have  introduced  a  micro-force  tensor  £  such  that  £  n  P  represents  the  power 

Jl  J  I  ±  ± 

density  expended  across  surfaces,  an  internal  micro-force  vector  7 r.  such  that  7 t  P  is  the 
power  density  expended  by  the  material  internally  (this  micro-force  accounts  for 
dissipation  in  the  material),  and  an  external  micro- force  vector  7.  such  that  7 P  is  the 
power  density  expended  by  external  sources.  After  the  introduction  of  the  Helmholtz 
free  energy  ip  =  ip^£..,D.,P.,P.^ ,  and  an  analysis  of  the  second  law  of  thermodynamics, 
the  following  constitutive  relationships  result. 
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'l  ~  dP 
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dP 
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(3.P 
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with  (3..  positive  definite  ,  where 


(10) 


The  form  of  the  Helmholtz  free  energy  that  is  chosen  to  mimic  the  behavior  of 
ferroelectric  single  crystals  is, 
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(11) 


The  specific  values  for  the  constants  that  are  used  to  model  barium  titanate  are  listed  in 
the  Appendix.  The  representation  of  this  free  energy  contains  a  few  important  material 
parameters  to  be  used  for  the  normalization  of  the  results  in  this  paper,  including  the 
spontaneous  polarization  PQ ,  the  spontaneous  strain  £() ,  the  critical  electric  field  for 

homogeneous  180°  switching  E  ,  the  characteristic  stress  ci.  =  EqPo  /  eQ ,  and  the 
domain  wall  length  scale  Z  . 

Numerical  Implementation 

Solutions  to  Equations  (l)-(ll)  can  be  obtained  numerically  using  the  finite  element 
method.  The  nodal  degrees  of  freedom  are  chosen  to  include  the  mechanical 
displacements,  the  polarization  components,  and  the  electric  potential,  from  which 
strains,  polarization  gradients,  and  electric  field  values  are  computed  within  the 
elements.  The  following  statement  of  the  principal  of  virtual  work  is  used  to  implement 
the  finite  element  method. 


fp.PSPdV  +  [  (a  ..Se  -  D  6E  +  r].6P  +  £JP  )dV 

J  V  3  »  J  \  Ji  V  i  %  'll  ^jl  l,j ) 
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(12) 


Standard  finite  element  procedures  lead  to  a  set  of  non-linear  algebraic  equations  for  the 
nodal  degrees  of  freedom,  and  the  solution  is  obtained  using  an  iterative  Newton- 
Raphson  method.  Additional  details  on  both  the  theoretical  foundations  and  the 
numerical  implementation  of  the  phase- field  model  can  be  found  in  References  [1-3]. 

Domain  Wall  Dislocation  Interactions 

Using  these  theoretical  foundations  and  numerical  methods  we  have  studied  the 
interactions  of  domain  walls  with  defects,  including  dislocations,  crack  tips,  and  free 
surfaces  and  interfaces  in  thin  films.  Here  we  present  results  from  applying  the  phase- 
field  framework  to  investigate  the  interactions  between  domain  walls  and  arrays  of 
dislocations  in  ferroelectric  single  crystals.  The  non-linear  finite  element  method  was 
used  to  determine  equilibrium  solutions  for  the  coupled  electromechanical  interactions 
between  a  domain  wall  and  a  dislocation  array.  The  numerical  simulations  demonstrate 
the  effect  of  the  relative  size  and  orientation  of  dislocations  on  180°  and  90°  domain  wall 
configurations.  In  addition,  results  for  the  pinning  strength  of  dislocations  in  the  case 
that  domain  walls  move  due  the  application  of  external  electric  field  and  shear  stress  are 
computed.  The  presented  numerical  results  were  compared  with  the  findings  reported 
for  charged  defects  and  it  is  shown  that  non-charged  defects,  such  as  dislocations,  can 


also  interact  strongly  with  domain  walls,  and  therefore  affect  the  ferroelectric  material 
behavior. 


Er/Ec=  0  Ev/Ec=  0.5  E/E=  1 


Fig.  1.  Polarization  distributions  in  the  y-direction  near  a  180°  domain  wall  interacting 
with  an  array  of  dislocations  in  the  [110]  direction.  The  x-scale  is  normalized  with  l()  and 

the  polarization  values  with  PQ .  The  solid  white  arrows  designate  the  direction  of  the 
polarization  vector.  Periodic  boundary  conditions  are  enforced  on  the  lines  at  half  the 
distance  h  =  100/()  between  the  dislocations.  For  E^/E  >1  the  wall  breaks  through 
the  dislocation. 

Figures  1  and  2  illustrate  the  types  of  results  produced  by  our  simulations.  Upon 
introducing  an  array  of  dislocations,  Fig.  1  demonstrates  the  effects  of  their  interactions 
with  the  180°  domain  wall  configuration.  Specifically,  P '  polarization  distributions  are 

plotted  near  dislocations  with  b  =  6[1 10]  and  b  =  10e()/f) .  The  distance  h  between  the 
dislocations  is  equal  to  100  Z() .  Fig.  1  shows  that  within  the  phase-field  framework,  the 

domain  wall  is  a  diffuse  zone  where  polarization  changes  smoothly  between  the  adjacent 
domain  states  which  are  marked  by  solid  white  arrows.  In  addition,  notice  that  the  wall 
is  kinked  with  respect  to  its  original  straight  configuration.  In  fact,  the  numerical 
results  reveal  that  as  the  dislocation  size  increases  the  kinking  becomes  more 
pronounced,  as  expected,  indicating  stronger  interactions  between  the  dislocations  and 
the  walls.  Domain  wall  kinking,  with  a  different  morphology,  has  also  been  reported  for 
domain  walls  interacting  with  arrays  of  point  charges.  These  results  demonstrate  that 
the  existence  of  non-charged  defects,  such  as  dislocations,  cause  similar  electroelastic 
interactions  in  ferroelectrics.  It  is  further  observed  in  Fig.  1  that  the  equilibrium 
position  of  the  domain  wall,  which  was  initially  placed  along  the  y- axis,  lies  to  the  left  of 
the  arrays  of  dislocations  when  no  electric  field  is  applied.  However,  as  the  applied 
electrical  field  increases  the  domain  wall  shifts  to  the  right  until  eventually  it  breaks 
through  the  array  of  dislocations.  In  Fig.  1,  E  is  the  critical  value  of  the  electrical  field 


required  to  force  the  domain  wall  to  break  though,  and  for  the  dislocation  size  shown  it 
is  equal  to  0.049F() .  It  is  interesting  to  observe  the  shape  of  the  kinking  and  note  that 

the  wall  always  bows  towards  the  dislocation  sites.  This  is  a  consequence  of  the  details 
of  the  electrostatic  fields  near  the  defects.  Hence,  it  is  apparent  that  dislocations  cause 
electromechanical  interactions  that  change  the  shape  of  a  180°  wall,  and  force  the  wall 
to  be  offset  from  the  array  when  no  external  loads  are  applied. 
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Fig.  2.  Polarization  distributions  in  the  y-direction  near  a  90°  domain  wall  interacting 
with  an  array  of  dislocations  in  the  [110]  direction  having  Burger’s  vectors  size 

b  =  lOe  The  x-scale  is  normalized  with  l  and  the  polarization  values  with  P  .  The 
solid  white  arrows  designate  the  direction  of  the  polarization  vector.  Periodic  boundary 
conditions  are  enforced  on  the  lines  at  half  the  distance  h  =  100Z()  between  the 

dislocations.  For  £  /  E  >1  the  wall  breaks  through  the  dislocation. 


Figure  2  shows  the  effect  of  applied  positive  electric  field  on  the  configuration  of  a  90° 
wall  interacting  with  an  array  of  dislocations  in  the  [110]  direction.  The  electric  field 

values  in  Fig.  2  are  normalized  with  the  critical  value  E  =  0.024 E  ,  which  causes  the 

wall  to  break  through  the  dislocations  for  a  Burger’s  vector  size  b  =  10e0/Q .  As  seen  in 

this  figure,  for  no  applied  electric  field  the  electromechanical  interactions  between  the 
domain  wall  and  the  dislocation  cause  significant  kinking  of  the  wall.  Compared  to  the 
180°  wall,  the  kinking  of  the  90°  wall  is  more  pronounced  for  the  same  dislocation  size. 
The  equilibrium  position  of  the  90°  wall  is  centered  the  array  of  dislocations  and  the 
application  of  positive  electrical  field  causes  the  domain  wall  to  shift  to  the  right.  When 
the  applied  field  reaches  the  critical  value  P  ,  the  results  shown  in  Fig.  2  illustrate  the 

equilibrium  position  just  before  the  wall  breaks  through  the  dislocation  array.  Similar 
results  have  been  obtained  for  a  90°  wall  interacting  with  [110] -type  dislocations.  It  is 
further  noted  that  the  effects  on  the  domain  wall  caused  by  its  interaction  with  the 


array  of  dislocations  in  the  case  that  electric  field  is  applied  are  very  similar  to  the 
corresponding  effects  caused  when  only  shear  stresses  are  applied.  This  obervat.ion 
agress  with  the  case  of  charged  defects  interacting  with  domain  walls.  Additional 
details  of  this  study  can  be  found  in  Reference  [2] . 

Domain  Nucleation  at  a  Crack  Tip 

Next,  results  for  the  growth  of  a  new  domain  from  a  crack  tip  during  purely  electrical 
loading  are  reported.  To  generate  the  solution  for  a  final  equilibrium  domain 
configuration,  the  domain  was  nucleated  at  the  crack  tip  and  then  allowed  to  evolve 
with  a  non-zero  polarization  viscosity  term  (i.e.  the  first  [5- term  in  Equation  (12)).  The 
loading  is  applied  by  first  ramping  up  a  uniform  charge  load  on  the  top  and  bottom 
surfaces  with  a  charging  rate  of  [3iju 4  /  E  =0.1  to  a  total  charge  of  u4  /  P  =  0.16  in  a 

60(.  x60(.  domain  with  an  electrically  impermeable  crack  half  of  the  way  through  its 
midplane.  The  charge  was  then  fixed  at  lua/Pq  =  0.16  and  the  domain  structure  was 

allowed  to  evolve  until  the  solution  was  sufficiently  close  to  the  equilibrium 
configuration,  at  which  point  the  polarization  viscosity  term  was  set  to  j3  =  0  to  find 
the  true  equilibrium  solution.  Thereafter  additional  charge  is  applied  to  the  surface  to  a 
final  value  of  c o ^  /  P  =  0.2  and  finally  the  charge  is  removed  from  the  surface  to  return 

to  the  initial  uncharged  state.  Note  that  in  order  to  ensure  accuracy  of  the 
computations  at  least  five  finite  element  nodes  span  any  domain  wall,  and  the  path- 
independence  of  the  ./-integral  is  verified  for  all  cases  of  equilibrium.  Details  of  the 
derivation  of  the  ./-integral  for  this  phase- field  theory  can  be  found  in  Reference  [3].  If 
the  mesh  is  too  coarse  then  mesh-pinning  of  the  domains  occurs  and  significant  but 
artificial  path-dependence  appears  at  equilibrium  in  the  /-integral.  Figure  3  shows 
contour  plots  of  the  y  component  of  the  polarization  distributions  at  four  different  times 
in  the  domain  evolution.  Figure  3  shows  that  the  90°  domain  needle  is  nucleated  at  the 
crack  tip  and  propagates  through  the  entire  domain  until  it  reaches  the  charged 
boundary. 

This  non-equilibrium  propagation  of  the  domain  supports  the  hypothesis  that  an 
instability  in  the  equilibrium  solution  exists  at  the  domain  nucleation  threshold. 
Additionally,  the  equilibrium  configurations  just  prior  to  the  domain  nucleation  with  no 
domain  (not  shown  in  the  figure)  and  the  domain  configuration  shown  in  Figure  3b  both 
occur  at  a  charge  loading  level  of  u>A  /  P  =  0.16,  and  are  sufficiently  distinct  from  one 

another.  This  unstable  growth  of  the  domain  is  in  contrast  to  domain  switching  zones 
predicted  using  phenomenological  constitutive  laws.  The  explanation  for  the 
discrepancy  is  that  these  phase-field  simulations  assume  a  defect-free  material.  In  such 
a  material  domain  walls  do  not  become  pinned  and  are  free  to  move  at  vanishingly 
small  levels  of  electromechanical  driving  force. 


Fig.  3.  Contour  plots  of  the  y  component  of  the  polarization  normalized  by  P0  in  the 
vicinity  of  the  crack  tip  for  (A)  ouA  /  PQ  =  0.16  during  the  non-equilibrium  evolution  of 

the  domain,  (B)  uj  4  /  PQ  =  0.16  at  the  final  equilibrium  state  for  the  domain,  (C) 
equilibrium  at  uj  A  /  P  =  0.2,  and  (D)  equilibrium  at  u) A  /  P  =  0.11 .  Only  the  upper 
half  of  the  region  is  displayed.  The  x  and  y  distances  are  normalized  distances  x  / l 
and  y  /  l  ,  and  the  polarization  scale  is  normalized  by  the  spontaneous  polarization  P  . 


Next,  results  for  the  apparent  crack  tip  energy  release  rate  calculation  are  presented. 
Figure  4a  plots  the  apparent  energy  release  rate  as  a  function  of  the  applied  charge 
loading  for  the  sample.  Note  that  points  A-D  on  Figure  4a  correspond  to  the  domain 
structures  illustrated  in  Figures  3A-3D  respectively.  Initially,  as  the  charge  is  applied 
the  energy  release  rate  is  negative  and  approximately  quadratic  in  the  applied  charge. 
These  features  of  the  energy  release  rate  are  in  accord  with  linear  piezoelectric  fracture 
mechanics  solutions.  Furthermore,  prior  to  the  nucleation  of  the  new  domain,  the 
solutions  are  carried  out  at  equilibrium  and  the  energy  release  rate  calculation  is  path- 
independent.  At  the  charge  load  level  of  u)  A  /  P  =  0.16,  nucleation  of  the  new  domain 
occurs,  and  at  this  point  the  crack  tip  energy  release  rate  is  approximately 
(g  =  — 3P(  P  /(j .  The  charge  is  held  fixed  at  this  point  and  as  the  new  domain  grows  the 

energy  release  rate  increases  from  G  =  —3 EPl  to  G  =  0.5 EPL  . 
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Fig.  4.  (a)  The  crack  tip  energy  release  rate  as  a  function  of  the  applied  charge.  Points 
A-D  correspond  to  the  domain  structures  illustrated  in  Figures  3A-3D  respectively,  (b) 
The  apparent  energy  release  rate  for  domain  structures  A  (blue,  non-equilibrium)  and  D 
(red,  equilibrium). 

Note  that  the  domain  evolution  at  loa  /  P  =  0.16  is  a  non-equilibrium  process  during 

which  the  apparent  energy  release  rate  calculation  is  path-dependent.  Specifically,  path- 
dependence  of  the  apparent  energy  release  rate  calculation  for  domain  structure  A  is 
plotted  as  the  blue  curve  in  Figure  4b.  Domain  structure  B  is  again  an  equilibrium 
configuration  and  the  energy  release  rate  calculation  is  path-independent.  After  domain 
structure  B  stabilizes,  additional  charge  is  applied  and  the  domain  structure  is  allowed 
to  evolve  at  equilibrium  to  domain  structure  C.  During  this  loading  process  the  energy 
release  rate  increases  in  an  approximately  linear  fashion.  Upon  reaching  structure  C  the 
applied  charge  is  removed  and  the  domain  structures  and  the  energy  release  rate 
“unloads”  along  its  original  loading  path  to  structure  B.  At  this  point  the  unloading 
path  diverges  from  the  original  loading  path  and  a  hysteresis  appears  in  the  energy 
release  rate  versus  applied  charge  response.  Domain  structure  D  is  arrived  at  during  the 
equilibrium  unloading  process  and  the  energy  release  rate  calculation  is  path- 
independent  as  shown  by  the  red  curve  in  Figure  4b.  Further  unloading  of  the  charge 
causes  the  domain  to  vanish  and  the  original  negative  quadratic  branch  of  the  energy 
release  rate  response  is  followed. 

The  most  interesting  aspect  of  this  simulation  is  the  departure  from  the  results  of  linear 
piezoelectric  fracture  mechanics.  Specifically,  this  calculation  is  the  first  that  we  are 
aware  of  that  predicts  that  the  crack  tip  energy  release  rate  can  be  positive  under  purely 
electrical  loading  for  impermeable  crack  face  boundary  conditions.  Furthermore,  the 
calculation  shows  that  an  existing  domain  structure  near  the  crack  tip  can  cause  a 
qualitatively  different  behavior  for  the  energy  release  rate,  positive  and  increasing  with 
applied  charge,  from  what  is  expected  in  linear  piezoelectricity,  negative  and  decreasing 
with  applied  charge.  We  do  note  that  large  scale  switching  does  occur  in  this  simulation 
and  so  a  direct  comparison  to  linear  piezoelectric  fracture  mechanics  concepts  is 


tenuous.  However,  these  simulations  demonstrate  the  effects  that  near  tip  domain 
structures  can  have  on  the  fracture  process  in  ferroelectric  crystals.  Specifically,  the 
negative  contribution  of  the  energy  release  rate  from  applied  electric  fields  may  in  fact 
be  positive  for  certain  domain  structures  near  crack  tips.  Hence,  the  modeling  of  crack 
tip  domain  structures  and  large  scale  domain  switching  behavior  in  fracture  specimens 
may  be  a  key  to  understanding  the  plethora  of  seemingly  disparate  experimental 
observations  on  electromechanical  fracture  of  ferroelectrics. 

Infinite  Boundary  Conditions  for  Crack  Problems 

During  our  studies  of  the  switching  zones  near  crack  tips  we  found  the  need  to  develop 
infinite  boundary  conditions  for  two-dimensional  crack  being  loaded  by  a  7b-field.  To 
accomplish  this  we  formulated  a  coupled  analytic/finite-element  method  for  two- 
dimensional  crack  problems.  With  the  method  two  classes  of  problems  can  be  studied. 
The  first  considers  problems  where  non-linear  constitutive  processes  occur  in  a  region 
near  the  crack  tip  and  the  remotely  applied  loading  can  be  characterized  by  the  linear 
elastic  ib-field  and  perhaps  the  T-s tress.  In  this  case,  the  finite-element  method  is 
applied  in  a  circular  region  around  the  crack  tip  where  non-linear  constitutive  response 
is  occurring,  and  stiffness  contributions  associated  with  a  numerically  implemented 
Dirichlet-to-Neumann  map  are  imposed  on  the  circular  boundary  to  account  for  the 
large  surrounding  elastic  domain  and  the  remote  applied  loading.  This  is  the  class  of 
problems  relevant  to  the  study  of  switching  zones.  The  second  class  of  problems 
considers  entirely  linear  elastic  domains  with  irregular  external  boundaries  and/or 
complex  applied  loadings.  Here,  the  discrete  Dirichlet-to-Neumann  map  is  used  to 
represent  a  circular  region  surrounding  the  crack  tip,  and  finite-elements  are  used  for 
the  external  region.  In  this  case  the  mixed  mode  stress  intensity  factors  and  the  T- 
stress  are  retrieved  from  the  map.  This  application  of  the  method  results  in  a  highly 
accurate  set  of  “crack-tip”  elements  that  can  be  used  to  decompose  the  stress  intensity 
factors  from  the  T-stress  in  linear  elastic  (or  linear  piezoelectric)  fracture  mechanics 
problems.  Full  details  of  the  method  can  be  found  in  Reference  [4], 

Path- Dependence  of  the  ./-Integral 

While  this  study  seems  somewhat  tangential  to  the  originally  proposed  work,  it  arose 
from  our  interests  in  how  the  ferroelast-ic  switching  zones  near  cracks  in  ferroelectrics 
affect  the  crack  tip  energy  release  rate.  The  section  just  prior  to  this  one  describes  our 
work  on  the  numerical  implementation  of  the  true  small-scale  yielding/switching 
boundary  conditions  for  such  crack  problems.  Our  first  test  of  the  method  was  simply 
to  look  at  the  path-dependence  of  the  J-integral  in  a  standard  elastic-plastic  material 
governed  by  J2-flow  plasticity  theory.  The  existing  body  of  decades  old  literature  on  the 
topic  suggested  that  the  J-integral  should  be  nearly  but  not  exactly  path-independent. 
What  we  found  in  this  study  was  interestingly  different  and  is  shown  in  Figure  5. 


Fig.  5.  Values  for  the  ./-integral  normalized  by  its  far-field  value  for  a  circular  contour 
of  radius  r  computed  by  the  domain  integral  method  near  a  crack  tip  under  mode  I 
loading  in  an  elastic-perfectly-plastic  material  with  v  =  0.3 .  The  markers  correspond  to 
different  points  along  the  load  history  and  different  sizes  of  the  plastic  zone  relative  to 
the  minimum  radial  dimension  of  the  elements  surrounding  the  crack  tip.  R  is  the 
approximate  plastic  zone  size  and  h1 .  is  the  minimum  element  size  near  the  crack  tip. 


The  results  shown  in  Figure  5  indicate  that  there  is  significant,  -20%,  path-dependence 
of  the  J-integral  as  the  radius  of  the  integration  path  is  reduced.  Additional  details  of 
this  study  can  be  found  in  Reference  [5].  From  an  informal  poll  we  determined  that 
several  experts  in  the  area  of  elastic-plastic  fracture  mechanics  found  this  results  to  be 
surprising.  Hence,  this  work  was  able  to  fill  a  gap  in  the  existing  literature  on  the  path- 
dependence  of  J  in  elastic-plastic  materials.  Continuing  work  in  progress  on  J  in 
ferroelast.ic  materials  has  shown  that  J  can  actually  be  elevated  at  the  crack  tip  due  to 
ferroelast.ic  switching. 
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Appendix 

Equation  (11)  presented  the  general  form  of  the  Helmholtz  free  energy  used  in  this 
paper.  For  a  coordinate  system  aligned  with  the  (100)  directions,  the  specific  form  used 
to  fit  the  dielectric,  piezoelectric  and  elastic  properties  of  ferroelectric  single  crystals 
that  undergo  a  cubic  to  tetragonal  phase  transformation  is  given  in  Equation  (Al). 
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In  the  equation  above  n0  =  8.854  x  10 


ox  =  -0.668325  EjPo 
Oi=12.4421£0/P0s 
^=  2.54138^/^ 
cx  =  2.04999  er(l/e0 
/,=  0-663581  P„/<P0 
4  =  0.687281 


Vm/C ,  and 

<*,  =  —3.80653  P0/P03 

a,  =  0.134226  A /p5 

6  0/0 

b2  =  1.74267  EjeaP0 
c2  =  0. 971673  <70/e0 
/2=  0.841326  P0/<P0 
/5  =  0.106647  PS/E(;P0 


oa  =  0.78922  B0/P/ 

6S  =  0.399353^/^ 

<=,  =  1-27976vA» 

/s=— 0.170635  B0/<P0 
/,=  0.213294  B0/<P„ 


s,  =  -3.66149  Eje0  P>  g ,  =  6.27423  F0>  9j  =  -1.21644^  i? 

where  er()  =  T/P  /  eQ  =  692x10°  N/m2.  In  addition,  PQ  =  0.26  C/m2 ,  eQ  =  0.0082  and 
PQ  =  2.182  xlO'  V/m  correspond  to  properties  of  monodonrain  single  crystal  barium 
t-itanate  at  room  temperature. 

The  parameter  a  appearing  in  Eq.  (Al)  determines  the  domain  wall  thickness.  If 
af)  =  1  x  10~10  V  •  nr  ’  /  C  a°  then  l  =  1  nnr  ,  and  therefore  the  180°  domain  wall  has 
thickness  equal  to  2  nnr  which  is  in  general  agreement  with  experimental  observations. 


